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Abstract 

We study the static and dynamic properties of bromine electrosorp- 
tion onto single-crystal silver (100) electrodes by Monte Carlo simula- 
tion. At room temperature the system displays a second-order phase 
transition between a low-coverage disordered phase at more negative 
electrode potentials and a c(2 x 2) ordered phase with bromine cover- 
age 1/2 at more positive potentials. We explore the phase diagram and 
demonstrate that the broad shoulder observed in room-temperature 
cyclic voltammograms is due to local fluctuations resembling ordered 
phases with coverage 1/4, which exist in the model at much lower 
temperatures. We construct a dynamic Monte Carlo algorithm using 
a thermally activated stochastic barrier-hopping model for the micro- 
scopic dynamics. We use this algorithm to study the phase ordering 
and disordering processes following sudden potential steps between the 
disordered phase and the c(2 x 2) phase, and to study the sweep-rate 
dependence in simulated cyclic- volt ammetry experiments. 

Keywords: Bromine adsorption. Continuous phase transition. Cyclic 
Voltammetry, Dynamic Monte Carlo simulation. Lattice-gas model. Potential- 
step experiments 
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1 Introduction 



A number of experimental techniques have emerged in recent years which 
enable studies of the structure and dynamics of electrode-electrolyte inter- 
faces in situ, i.e., without removing the electrode from the liquid IQ. This 
new atomic scale information allows greater understanding of the deposi- 
tion process where, in contrast to high-vacuum surface-science methods, the 
electrochemical environment provides direct experimental control over the 
electrochemical potentials of the adsorbates, rather than over their coverages 
or fluxes. These new in situ experimental techniques, together with new 
methods to manipulate adsorbates at the interface on the nanometer scale 
[|2|, promise to make electrodeposition an even more important method for 
low-cost, high-quality materials synthesis than it is today These develop- 
ments emphasize the importance of developing theoretical and computational 
methods to study the dynamical aspects of electrodeposition in microscopic 
detail. 

Motivated by these developments in electrochemical surface science, we 
here present a static and dynamic Monte Carlo (MC) simulation study of 
a simple system which can serve as a prototype for more complicated ones 
of applied interest. This system is a submonolayer of Br electrodeposited 
onto a single-crystal Ag(lOO) electrode. Recently it was investigated at room 
temperature by in-situ surface X-ray scattering (SXS) and chronocoulometry 
by Ocko, Wang, and Wandlowski That study detected a second-order 

phase transition between a disordered low-coverage phase at more negative 
electrode potentials and a doubly degenerate c(2 x 2) ordered phase with 
coverage 1/2 at more positive potentials. (This phase has also been observed 
by low-energy electron scattering (LEED) for Br adsorbed on Ag(lOO) in 
vacuum The transition was shown to belong to the two-dimensional 

Ising universality class [Q. The scattering data 0, Q ^-Iso indicate that Br 
adsorbs at the four-fold hollow sites between the surface Ag atoms, consistent 
with recent in situ X-ray absorption fine structure (XAFS) results [0. The 
fluctuations around the adsorption sites are small, suggesting that a lattice- 
gas approximation for the adlayer structure and dynamics is appropriate. 
Here we investigate in detail such a lattice-gas model with a nearest-neighbor 
excluded-volume interaction and lateral dipole-dipole repulsive interactions 
IP). The dipole interaction strength and the electrosorption valency, which 
are parameters of the model, are estimated by fitting simulated equilibrium 
adsorption isotherms to experimental chronocoulometric data 



2 



Cyclic voltammograms (CVs) at room temperature show two main fea- 
tures 0, 1^, 1^, |10|, 0: a broad shoulder at more negative potentials which 
ends in a sharp peak at more positive potentials. By SXS it was confirmed 
that the sharp peak corresponds to a second-order phase transition in the 
layer of adsorbed Br 0]. However, the cause of the broad shoulder has been 



less clear. It was previously suggested [|10] that interactions with surface 
water are responsible for the shoulder, but equilibrium MC simulations by 
Koper p, |12|, as well as our own work, show the same features in simulations 
without water. We show here that the broad shoulder in the simulated CV is 
due to fluctuations within the disordered adlayer. Locally, these fluctuations 
resemble other ordered phases with coverage 1/4, which exist in the model 
at much lower temperatures. 

Going beyond these equilibrium simulations, we further present a dynamic 
MC study of the ordering and disordering processes between the disordered 
phase and the c(2 x 2) phase, assuming a perfectly stirred solution so that 
bulk diffusion effects can be ignored. For ordering, we expect similar behavior 
to the random sequential adsorption process of the hard-square model with 
lateral diffusion |]13|, |14|. At late stages of the ordering, this dynamic model 
is expected to be in the same dynamic universality class as the ferromagnetic 
15| or antiferromagnetic Ising model in zero magnetic field following 



sudden cooling from infinite temperature to a temperature below the critical 
point. These phase ordering dynamics with non-conserved order parameter, 
often called Model A in the literature on dynamics of phase ordering and 



phase separation [0, |T^], are driven by the reduction of interface curvature 



[|18| , |19[. We further present dynamic results for the disordering process 
and discuss how simulation results for this process can be compared with 
experiment. 

We also apply the dynamic MC algorithm to study sweep-rate dependence 
in simulated CV experiments. Results are presented for the response of the 
system to a wide range of sweep rates, and we suggest some comparisons 
between simulation results and experiment. Some preliminary results of this 
study were included in Refs. pT |. 

The remainder of this paper is organized as follows. In Sec. |^ we introduce 
the lattice-gas model and the relation between the electrochemical potential 
and the experimentally measured electrolyte concentration and electrode po- 
tential. In Sec. 1^ we present our equilibrium results, and in Sec. ^ we intro- 
duce our dynamic MC algorithm. In Sec. ^ we present dynamic simulations 
of sudden potential-step experiments, and in Sec. ^jwe present dynamic simu- 
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lations of CV experiments. In Sec. |^ we give a summary and our conclusions. 
The Appendix describes the nonhnear fitting procedure used to determine 
the lattice-gas parameters. Animations and additional figures are available 
on the World Wide Web [||. 



2 Lattice-Gas Model 

The lattice-gas model consists of a two-dimensional L x L square array of Br 
adsorption sites, corresponding to the four-fold hollow sites on the Ag(lOO) 
surface. Periodic boundary conditions were used to eliminate edge effects in 
the relatively small systems it is possible to simulate. The grand-canonical 
lattice-gas Hamiltonian is 



H = - (pijCiCj - fi^Ci , (1) 

i<j 1=1 

where i and j denote sites on the lattice, q denotes the occupation state of 
site i, which is either (empty) or 1 (occupied), Si<j denotes a sum over all 
pairs of sites on the lattice, (pij denotes the lateral interaction energy of the 
pair and JI is the electrochemical potential. The sign convention used 

here is such that 4>ij < denotes a repulsive interaction, and JI > favors 
adsorption. We measure and fl in units of meV/pair and meV/particle, 
respectively. [For brevity we write both simply as meV.] 

The electrochemical potential Jt is the intensive variable which controls 
the average number of Br in the adlayer. In the weak-solution approximation 
it is related to the electrode potential by 

Ti = -po + kBT\nj^-ejE , (2) 

where /Iq is an arbitrary reference level, k-^ is Boltzmann's constant, T is the 
absolute temperature, e is the elementary charge unit, [C] is the concentra- 
tion of Br~ in solution, [Cq] is an arbitrary reference concentration (here we 
use Co = 1 mM), 7 is the electrosorption valency (76 is the charge trans- 
fered through the external circuit for each Br ion adsorbed on the electrode 



23| , ^% P5|), and E is the electrode potential in mV. 

A variety of lateral interactions have been previously explored for this 
model, including finite and infinite nearest-neighbor repulsion and screened 
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dipole-dipole interactions From these previous studies and our own ex- 
ploratory simulations, we find that the interactions are quite adequately de- 
scribed with a nearest-neighbor excluded- volume interaction (infinite nearest- 
neighbor repulsion) related to the large ionic radius of Br" (1.96 A l2|]), plus 
a long-range dipole-dipole repulsion. The pairwise lateral interactions used 
here are 

{—oo r = 1 
V2<r<5 , (3) 
r > 5 

where r is the separation of an interacting Br pair within the adlayer, mea- 
sured in units of the Ag(lOO) lattice spacing (a = 2.889 A and 0nnn 
(which is negative) is the lateral dipole-dipole repulsion between next-nearest 
neighbors. In other work it has been noticed that reducing the infinite 
nearest-neighbor repulsion to finite but large values has very little effect on 
the isotherms All interactions are truncated beyond five lattice spacings 
to simplify computation. This truncation can be seen as an approximate 
cut-off due to screening. With this cut-off, the total interaction energy of a 
fully ordered c(2 x 2) phase is approximately 87% of the interaction energy 
without a cut-off. 

The Ag(lOO) substrate provides a square adsorption lattice of hollow sites 
for the commensurate Br adlayer. The nearest-neighbor excluded- volume in- 
teractions give rise to two different sublattices which have a geometry corre- 
sponding to the two colors on a checkerboard. Each of the two degenerate 
c(2 X 2) phases can occupy only one of the two sublattices, which we label 
A and B. The A sublattice coverage is the fraction of occupied sites on sub- 
lattice A, defined as 9^ = Nj^ J2ieA^i^ where Na is the number of sites on 
sublattice A and Z^ieA denotes a sum over all sites on that sublattice. The 
B sublattice coverage, 0^, is analogously defined. For this system, there 
are two linearly independent observables of interest: the total Br coverage, 
= (©A + 0_b)/2, and the staggered coverage, 6s = Qa — ©b, which is the 
order parameter for the c(2 x 2) phase. 

3 Equilibrium Monte Carlo Simulations 

Monte Carlo simulations were performed at room temperature {k-QT = 25 meV, 
corresponding to T ^ 290 K) to determine the equilibrium properties of the 
Br adlayer for different values of JI. For these equilibrium simulations, we 



5 



used the Metropolis algorithm with single-site updates at randomly se- 
lected sites. 

The main purpose of these equilibrium simulations was to obtain adsorp- 
tion isotherms for parameter estimation by fitting to experimental data. The 
most important region of comparison is in the disordered phase, where the 
fluctuations are small. It was therefore sufficient to use relatively small sys- 
tems with L =32. No significant differences were noticed between data for 
L = 32 and L = 128 (see below). No systematic finite-size scaling was per- 
formed near the phase transition. Thus, the results obtained in the critical 
region are only approximate. 

For a given temperature and chemical potential, we generated a set of 
samples from which we computed the averages (9(T, /!)) and (|9s(T, 7l)|), 
where (■) denotes an average over a set of statistically independent equilib- 
rium samples. For convenience, we drop the (■) notation for the equilibrium 
isotherms. 

The critical value of the electrochemical potential and the shape of the ad- 
sorption isotherm in the disordered phase are quite sensitive to the strength 
of the long-range repulsive interactions. To determine the dipole interac- 
tion strength, 0nnn, and the electrosorption valency, 7, we therefore simul- 
taneously compared the simulated room-temperature adsorption isotherm, 
0(/2), with experimental isotherms, Q'^^'^{E), obtained from chronocoulom- 
etry ||^. The nonlinear fitting method is described in the Appendix. The 
resulting parameter values are (^nnn = —26 ± 2 meV, 7 = —0.73 ± 0.03, 
and fiQ = 0.67 ± 0.03 meV. These results are consistent with those found 
previously by others [|, |, |] 

The fitted MC isotherm is shown together with the chronocoulometric 
coverage isotherms in Fig. |l|. For the two lowest concentrations, the agree- 
ment is very good over the whole range of electrode potentials, while the 
agreement in the ordered phase is somewhat less satisfactory for the highest 
concentration. This disagreement at higher concentrations suggests a break- 
down in the weak-solution approximation and the ultimate need for a more 
detailed treatment of the solution. 

Figure || shows the simulated Q{fi) and |9s(/i)| at T ^ 290 K, with 
0nnn = "26 mcV as obtained from the fits to the experimental isotherms. 
The order parameter, Bs, is zero in the disordered phase, while it is near ±1 
deep in either of the two degenerate ordered phases. A second-order phase 
transition between the low-coverage disordered phase and the c(2 x 2) ordered 
phase is observed at a critical electrochemical potential, flc ~ 180 ± 5 meV. 
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As p, |0s| oc {fx — ficY^^ 0], as expected for the Ising universality 

class. 

The coverage isotherm shows a long, gradual slope in the disordered 
phase, which is not present in a model with only nearest-neighbor exclusion 
(called the hard-square (HS) model) [^, As a result, fic is considerably 
larger than that of the pure hard-square model, ~ 33 meV ||2^, 29[]P]. 
The simulated and experimental isotherms show coverages at ftc of about 
0.37 ± 0.01 and 0.36 ± 0.02, respectively. Within statistical errors, this is 
the same result as the pure hard-square model, 6^^ = 0.368 ^30- The 
adsorption isotherm, 0(/i), has its maximum slope at flc- Corresponding 
maxima in the slopes are also seen in the experimental isotherms in Fig. |l|. 
Simulation data are included in Fig. |^ for both L = 32 and L = 128. On the 
scale of this figure, differences between the two sizes are visible only in the 
data for |6s| very close to fic, where d|9s|/d/2 is expected to diverge with 
increasing L as L^/"^ ||31|| P|. 



A quasi-equilibrium CV (the CV in the limit of vanishing potential sweep 
rate) can be obtained from the simulated equilibrium isotherm using the 
simple relation, 

,de , ded/idE 2 2,.dedE de 

where JT" is the voltammetric current density (oxidation currents positive), 
M is the total number of adsorption sites per unit area (a~^ = 1.198 x 10^^ 
sites/cm^ 0]), and dE/dt is the sweep rate. Figure | shows d0/d/2, which 
exhibits the same broad shoulder and sharp peak seen in experimental CVs 
IQ, 1^, 0. The shoulder has previously been attributed to interactions with 



surface water |T0[ . However, the fact that the lattice-gas model does not con- 
tain explicit water molecules indicates that the shoulder is most likely caused 
by configurational fluctuations in the Br adlayer, as was also pointed out by 
Koper p[. Differences between the simulation data for L = 32 and 128 are 



^ The value fL^^ ~ 33 meV is obtained as p,^^ = ksTlnZc, using precise numerical val- 
ues for the critical activity of the hard-square model, Zc = 3.7962(1) and Zc — 3.7959(1) 
[p9| , obtained by series expansions and transfer-matrix finite-size scaling, respectively. 

^ Precise numerical results have been obtained for the pure hard-square model by 
transfer-matrix calculations, including 9^ = 0.36776(1) [|§ and 6^ = 0.368(1) ||. 

^ At p,c in a finite system of size L, d|8s|/d/2 cx L''^'^, where 7 and v are the critical 
exponents for the ordering susceptibility and the correlation length, respectively. For the 
two-dimensional Ising model 7 = 7/4 and 1/ = 1 p^, |li. 
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only visible near the peak, which is expected to diverge logarithmically with 
L [^0. This quasi-equilibrium CV is compared with dynamically simulated 
CVs for different potential sweep rates in Sec. ^. 

Having obtained an estimate for the lateral interaction constant by equi- 
librium simulations at room temperature, we briefly explore the full phase 
diagram of the model with this interaction strength, shown in Fig. ^. We 
expect the phase diagram for this system to be topologically similar to that 
for the hard-square model with next-nearest neighbor repulsion |Q (even 
though the low-coverage ordered phases in the two systems are different). In 
addition to the disordered phase and the two equilibrium ordered phases [the 
low-temperature c(4 x 2) phase with B = 1/4 and the c(2 x 2) phase with 
= 1/2 which is observed at all temperatures for sufficiently large /i], there 
is also a metastable phase, p(2 x 2) with B = 1/4, which lies closely above 
the c(4 X 2) phase in energy. To facilitate equilibration at low temperatures. 



in this simulation we used a heat-bath algorithm which updates clusters |34 



consisting of randomly chosen nearest-neighbor pairs of sites. The presence 
of the p(2 X 2) phase made low-temperature equilibration very difficult, even 
with the pair- update algorithm. Other phases with < 1/4 probably also 
exist as a consequence of the long-range repulsion; however, we expect these 
phases to be present only at extremely low temperatures, and their effects 
should be negligible at room temperature. The phase boundaries have been 
determined by ground-state calculations for T = 0, by the location of the 
sharp peak in the quasi-equilibrium CV for T > 100 K, and by inspection of 
the pair-update MC sample configurations for < T < 100 K. 

The broad shoulder in the CV is caused by configurational fluctuations in 
the disordered phase, which locally resemble the 6 = 1/4 phases [^, At 
the experimentally unrealistic low temperatures where the c(4 x 2) phase is 
observed in the model, the shoulder breaks up into two transitions between 
the disordered phase and the c(4 x 2) phase. Although the c(4 x 2) phase 
is not experimentally observable in this system, local fluctuations which re- 
semble the = 1/4 phases can be seen in the model in equilibrium at room 
temperature and should be observable in the diffuse intensity of SXS exper- 
iments. We therefore next present simulated equilibrium SXS intensities in 
the room-temperature disordered phase at B ~ 1/4. 

At flc d8/d/i is expected to be proportional to the heat capacity of the adlayer. 
The heat capacity of the two-dimensional Ising model diverges logarithmically with L 
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We define the structure factor as S{k) = |c(fc)p, where |c(fc)p denotes 
the squared absolute value of the complex Fourier transform of the real-space 
adlayer configuration {q}, and k = [k^, ky) is the wave vector in the plane of 
the adlayer (given in units of the inverse lattice constant, a~^). Equivalently, 
S{k) can be expressed as the Fourier transform of the adsorbate correlation 
function, C{fi — fj) = (ciCj). The Fourier transform is normalized such that 
S{k = (0,0)) = (LG)^. In the Born approximation, the structure factor 
is proportional to the SXS intensity In Fig. ^(a) and (b) we show 

the simulated structure factor for the disordered room-temperature phase at 
/i = +100 meV, where 6^1/4. Small domains of c(2 x 2), c(4 x 2), and 
p(2 X 2) symmetry are noticeable in snapshots of the adsorbate configurations 
[Fig. ^(c)], and they also affect the simulated diffuse scattering intensity. 
Referring to Fig. |^(b), Bragg peaks would occur at the points marked x 
for a fully ordered c(2 x 2) phase, at the points marked Q fo^ P(2 x 2), 
and at the points marked □ for c(4 x 2). Due to the disorder, all of these 
peaks [except for the one at A; = (0, 0)] are very broad and combine to 
give the anisotropic diffuse scattering intensity shown in the figure. The 
individual peaks would become narrower and more distinct at temperatures 
less than room temperature (but still above the critical temperature for the 
low-temperature ordered phase) or with stronger dipole-dipole interactions. 



4 Dynamic Monte Carlo Simulations 

Next we investigate the dynamics of the Br adlayer using a dynamic MC 
algorithm. The microscopic processes considered are adsorption/desorption, 
in which a particle enters or leaves the adlayer, as well as nearest-neighbor and 
next-nearest neighbor lateral diffusion. The thermally activated, stochastic 
barrier-hopping picture of the microscopic dynamics used here is described 
in further detail elsewhere [0S . 



We assume a simple linear form for the local energy landscape, which 
is equivalent to the Butler- Volmer approximation with symmetry parameter 
1/2 [^ |39|. See Fig. |. The initial and final states, / and F, respectively, 
are represented by lattice-gas configurations and differ by a single micro- 
scopic move. We label each type of microscopic move (adsorption, desorption, 
nearest-neighbor and next-nearest neighbor diffusion) by the index A. The 
state T\ represents an intermediate transition state of higher energy between 
/ and F, which is associated with the move A. This intermediate state can- 
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not be represented by a lattice-gas configuration, and we associate with it a 
"bare" barrier A^. The energy of the transition state is then approximated as 



Ut^ = + Ax , (5) 

where Ut^ , Uj, and Up are the energies of the states Tx, I, and F, respectively. 
Although other choices of the transition probability are used in the liter- 



ature as well |40], we here approximate the probability per MC step per site 



(MCSS) of making a transition from / to F by the simple one-step Arrhenius 



rate m, O, 42 



^ .exp (-^) ^ .e.p (-^) exp {-"-L^) . (6) 

where the rate constant u (in units of (MCSS)^^) sets the overall time scale 
for the simulation and is assumed to be the same for all microscopic processes. 

Dynamic Monte Carlo steps proceed in the following way. First, a lattice 
site j is randomly selected. If j is empty, only adsorption may be attempted. 
If j is occupied, both desorption and lateral diffusion may be attempted. 
In general, some moves will be excluded by the presence of other adatoms 
(corresponding to Up = +00). Next, a weighted list is constructed of new po- 
tential configurations F which would result from successfully attempting each 
of the allowed moves for site j. The probability to undergo each transition 
is then calculated from Eq. (^. Blocked moves have transition probability 
zero. To ensure microscopic reversibility (detailed balance), the possibility 
that no change in configuration occurs {F = I) must also be included in the 
list of new configurations. The probability to remain in the initial state is 
R{I\I) = 1-Ef^jR{F\I). 

In general, there is a possibility that R{I, I) could be negative. However, 
this essentially means that the time corresponding to one MCSS is chosen 
too long for a consistent stochastic description. Thus, z/ must be chosen 
sufficiently small that -R(/|/) remains positive. As long as this condition is 
satisfied, changing u only corresponds to changing the overall time scale for 
the simulation. In our simulations we used u = 1 (MCSS)~^, which ensured 
that R{I\I) remained positive with the barriers used here. 

The main price paid for the simplicity of coding our algorithm is that it is 
rather slow in parameter regimes where the acceptance probabilities are very 
low. However, except for extremely fast processes, in which the acceptance 
probabilities are large and change significantly within a single MCSS, it is 
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1 0, 



equivalent to more sophisticated continuum-time algorithms |¥T, 

mm- 

The adsorption/desorption barrier is Aa, the nearest-neighbor diffusion 
barrier is Ann, and the next-nearest neighbor diffusion barrier is Annn- In the 
simulations, the diffusion barriers are Ann — 100 meV and Annn — 

200 meV, 

based on ab-initio calculations |58] of binding energies for a single Br ion 
to a Ag(lOO) substrate in vacuum. The barrier Ann is approximated by the 
difference in binding energy between the four-fold hollow site and the bridge 
site, while Annn is approximated by the difference in binding energy between 
the four-fold hollow site and the on-top site. 

Theoretical estimates of Aa are extremely sensitive to the details of both 
the ion-surface and ion- water interactions Based on data from cal- 

culations of the potential of mean force for halide ions in water near a 
Cu(lOO) surface our very rough estimates indicate that values between 
200 and 500 meV are not unreasonable. Since we also expect Aa to be 
significantly larger than Annn, we examined some aspects of the behavior 
for Aa = 300 meV and Aa = 400 meV. Finding few qualitative differences 
between simulation results for these two values, we chose the lower value, 
Aa = 300 meV, to increase the simulation speed. 



5 Simulated Potential Steps 

One way to probe the dynamics of an adsorbed layer is to subject it to 
a sudden change in electrochemical potential. Dynamic simulations of such 
potential-step experiments are started by equilibrating the room-temperature 
system at a potential fii, using the simulation technique described in Sec. 3. 
Then, at time t = 0, the potential is instantaneously stepped to a new value 
/22, and the simulation continues by the dynamic MC algorithm described 
in Sec. ^. To study the dynamics of the phase ordering and disordering 
processes, we examined four different potential steps: two from the disor- 
dered phase into the ordered phase, and two from the ordered phase into the 
disordered phase. 

Before discussing each of these potential steps in detail, we give the ini- 
tial adsorption, desorption, and diffusion rates at fi2, immediately after the 
sudden potential change. We define the probability during one MCSS of 
accepting (when possible) an adsorption move as i?a, a desorption move as 
i?d, a nearest-neighbor diffusion move as i?nn, and a next-nearest neighbor 
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diffusion move as -Rnnn- As an example, we present the rate calculations for 
a deep disorder-to-order step with fti = —200 meV and ^2 = +600 meV. 
Each of the four rates is calculated from Eq. (|^). The initial and final en- 
ergies, Uj and Up, respectively, are calculated from Eq. (p. At 9 ~ 0; 
therefore, at early times the surface is nearly empty. For adsorption, / is 
approximated as having 6 = 0, and F contains only one occupied site. 
Thus, Ui = 0, Up = -/i2, and = uexp [(-A^ - {Up - Ui)/2)/k^T] = 
uexp [0.0]. For desorption, I must contain at least one occupied site for the 
process to be possible; therefore, / is approximated as having only one oc- 
cupied site, and F is approximated as an empty surface. Thus, Uj = —fi2, 
Up = 0, and = z/exp [(-A^ - (Up - Ui)/2)/k^T] = z/exp [-24.0]. For 
nearest-neighbor diffusion, I is also approximated as having only one oc- 
cupied site, and F is one of the four possible configurations which would 
occur from a nearest-neighbor diffusion move. Thus, Uj = Up = —fl2, and 
Rnn = 4i/exp [(-Ann " (Up - Ui)/2)/kBT] = 4z/exp [-4.0], where the fac- 
tor of 4 comes from adding the contributions from each of the four final 
states. Similarly, for next-nearest neighbor diffusion, Ui = Up = —fi2, and 
Rnnn = 4:U Bxp [{- A^nn - {U p - U i) / 2) / ksT] = 4i/ exp [-8.0] . The initial 
rates for all four potential steps are summarized in Table 1. Note that the 
rates change significantly during the course of the time evolution following 
each step. 

We now discuss each of the four potential steps. First we consider a 
disorder-to-order step with fii = —200 meV and fl2 = +600 meV. At 
0, and the final potential fi2 is deep into the ordered phase, approx- 
imately 420 meV on the positive side of the transition at ftc ~ 180 meV. 
Figure 0(a) shows time series for both O and |9s|. For deep steps like these, 
the desorption rate is negligible (see Table 1). Thus, the adsorption dynam- 
ics are essentially described by the process of random sequential adsorption 
with diffusion (RSAD) |14|. The coverage quickly reaches the jamming 
coverage for random sequential adsorption of hard squares without diffusion, 
0j ~ 0.364 |5y, |51], This is only slightly less than the critical cover- 



age, 6c ~ 0.368, which is reached in approximately 37 MCSS. During this 
time, small, uncorrelated domains of the ordered phase are formed on both 
sublattices. As the lattice fills up, the interfaces between these domains con- 
tain empty sites, but the nearest-neighbor exclusion generally prevents the 



^ Precise values for Qj in the hard-square model, 0.364 13(1) and 0.364133 0(5), have 
been obtained by MC simulations |pfl and series expansions ||5]|, respectively. 
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adsorption of additional Br ions. This is followed by the growth of adsorbate 
domains on the two sublattices which is well described by dynamical scaling 
(see below). At this stage, the coverage can increase only where domain 
walls move together, opening a gap large enough to fit an additional Br. In 
the absence of single-particle diffusion, the domain walls can move only by 
a series of very improbable desorption/adsorption events. However, nearest- 
neighbor diffusion allows some Br ions to diffuse between adjacent domains, 
thus providing a much faster mechanism for interface motion. For deep steps 
with the barriers used here, interface motion is almost exclusively due to 
nearest-neighbor diffusion. Excluding very small fluctuations on short time 
scales, the coverage is always driven toward the equilibrium value, 6eq = 1/2. 
Thus, the total length of interfaces is steadily decreasing. Finally, selection 
of one or the other of the two degenerate sublattices takes place at late times. 
The phase-selection time is expected to depend on L. 

In the dynamical scaling regime, the correlation length ^ of the system 



should grow as a power law with time t |15 



i oc t\ (7) 

In models with non-conserved order parameter such as the one studied here, 
the domain growth is driven solely by the tendency to reduce the amount 
of interface. For such systems (known as Model A in the nomenclature of 
critical dynamics [|r^, [T^) n = 1/2 |T^, |T^, |TB], 0. We have used two 



different measures for ^. One common measure uses the width of the diffuse 
structure-factor peak corresponding to the ordered phase, here the (1/2, 1/2) 
OT k = (tt, vr) peak. We define the correlation length given by the structure 
factor as ^5 = 27r/ (k), where 

(k) = ' "'"'^ . (8) 

Here {S{\k\)) is the circularly averaged S{k) and Z]|^_^ |>o ^^'^^ 
lowed values of \k\ ^ 1^,^,77 1, such that — /Ct^.ttI is less than n to avoid picking 
up contributions from the small diffuse scattering surrounding the k = (0,0) 
Bragg peak. Figure § shows three-dimensional views of the structure factor 
and typical system configurations at three different times following this deep 
step into the c(2 x 2) ordered phase. The diagonal anisotropy in S{k) at late 
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times corresponds to domain walls oriented at 45° to the lattice directions, 
which are stable under diffusion. Figure ^ shows the circularly averaged 
structure factor at various times. At late times the structure factor is pro- 
portional to I A; — /c7r,7r|~^. This asymptotic behavior, known as Porod's law 
[]53| , |54| , indicates that the domain walls are microscopically thin. 

For this system, we are also able to exploit the area covered by the empty 
domain walls as a second measure of C, [|T3l , , since a decrease in interfacial 
length is equivalent to an increase in the correlation length ||55[. We define 



this second measure as = (©eq — where Geq is the equilibrium value 
of the coverage at /X2 (here, 9eq = 1/2). The inset in Fig. |^(a) shows ^5 and 
as functions of t^/^ for L = 256. The agreement with the expected t^/^ 
behavior is excellent. This behavior has also been found in previous studies 



of RSAD ||T^, |1J] . For an infinite system, ^ would continue to grow without 
bound; however, when ^ reaches the same order of magnitude as L, the 
system enters the phase-selection regime in which one of the two degenerate 
ordered phases is randomly selected. 

For the other disorder-to-order step simulation, we take /2i = —200 meV 
and /i2 = -1-250 meV. This step ends approximately 70 meV on the positive 
side of the transition. Figure ^b) shows the dynamic behavior following 
this shallow step into the ordered phase. As one would expect, the overall 
response is much slower than for the deep step. Although the behavior at 
early times is more complicated than for the deep step, ^5 does show t^/^ 
domain growth at late times (see inset). The length scale evaluated from the 
coverage, ^w, is not an accurate measure for shallow steps in the time range 
we have simulated, since the interfaces are not as clearly defined as in the 
deep step. 

As expected, the ordering dynamics are well described by the dynami- 
cal scaling prediction of Model A. Although dynamical scaling has yet to 
be experimentally measured in this system, it has been measured in other 
electrochemical adsorption systems p6[, and we believe that t^^^ dynamical 



scaling will be experimentally observed in the future for this system as well. 

Potential steps that induce disordering are also of interest. For deep 
negative steps, fii=+600 meV and /i2=— 200 meV, the disordering process 
is relatively uneventful, as shown in Fig. pi]|(a). Particles simply desorb at 
a roughly constant rate, leading to an essentially exponential relaxation to 
equilibrium. 

Much more interesting behavior is seen after a shallow negative step, such 
as /2i=-|-600 meV and /22=-|-100 meV, shown in Fig. [TO|(b). The simulation 
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begins with all particles on sublattice A and relaxes to the disordered phase 
at j22, where O 1/4. Four different dynamical regimes are seen. In the first 
regime, particles simply desorb from the surface, so that dBs /dt ~ 2 dO / dt. 
In the second regime, both desorption and diffusion contribute significantly 
to the disordering process. As more sites become vacant by desorption, small 
domains of the B phase are created, both by lateral diffusion from sublattice 
A and by re-adsorption from the solution. As a result, the diffuse scattering 
peak narrows. In the third regime, the net desorption slows down, and 
diffusion and to some extent adsorption/desorption are the only significant 
contributions to the disordering process. Here, dGs/dt is largest. In the 
final regime, the two sublattices become approximately equally populated, 
yielding 0s 0. Similar behavior is seen in a dynamic mean- field theory 



The inset in Fig. |TD|(b) shows the correlation length obtained from the dif- 
fuse (1/2,1/2) scattering peak as a function of time for the step to +100 meV. 
This length can be interpreted as the typical distance over which the order- 
parameter fluctuations are correlated. Thus, at early times, when the Br 
are predominantly on sublattice A, the fluctuations onto sublattice B are 
very small. We believe that the temporary increase in the correlation length 
obtained from the diffuse scattering is caused by the system passing near 
the critical point of the phase transition. The length scale of the fluctua- 
tions decreases again as the coverage falls below its critical value and the 
macroscopic order is lost. 

For shallow disordering steps, time series for B and 0s are quite sensitive 
to the choice of the adsorption/desorption barrier. Here we discuss how this 
dependence can be exploited to estimate both Aa and the rate constant u 
from future experimental results. For simplicity, we assume that our choices 
of diffusion barriers are realistic. Neglecting this assumption does not affect 
our qualitative conclusions since the time scale for the disordering process 
is dominated by Aa — Ann- We first discuss 6(t'), which depends on time 
only through the dimensionless rescaled time t' = tz/exp (— Aa/A^B?"). See 
Fig. [TT|. By comparing experimental measurements of G(t) with the rescaled 
simulation results, the product z/exp (— Aa/ZcBT") can be estimated. However, 
knowledge of B(t) alone is not enough to independently estimate both u and 
Aa. To estimate both parameters, knowledge of |0s(^)| is also required. The 
dependence of |9s| on Aa can be seen in Fig. |lT], but the dependence is 
more readily seen in a parametric plot of |0s(0)|, which is independent of 



u. Such a plot is shown in Fig. 12. In this form, the shape dependence on 
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Aa is more obvious. Also, the result that d|0s|/dt = 2dQ/dt at short times 
is easily seen in the parametric plot. By plotting |Bs(G)| from experiments 
and simulations together, it should be possible to estimate Aa, assuming 
that the diffusion barriers used are realistic. Once Aa has been estimated, 
the product uexp (— Aa/fceT), known from comparing the time scales of the 
experimental and simulated processes, can be used to estimate u. It is our 
hope that experimental measurement of both Q{t) and |6s(^)| will be possible 
in the future for this shallow disordering process. 



6 Simulated Cyclic Voltammograms 

Next we use the dynamic Monte Carlo method described in Sec. ^ to study 
the dynamic response of the adlayer to a continuously varying potential, as 
realized in a CV experiment. The system displays a phase lag between the 
ordering and disordering processes and the time dependent jl. For these sim- 
ulated CVs, fL is linearly ramped from fii to fi2 with a sweep rate p and then 
ramped back to fii at the same rate. A wide range of sweep rates were simu- 
lated, from 10^^ to 3 X 10"^ meV/MCSS. Figure |T3| shows 6 as a function of fl 
and is analogous to the familiar hysteresis plots for magnetic systems, show- 
ing the magnetization as a function of the applied magnetic field [^. The 
net flux of Br~ into the adlayer is proportional to the voltammetric current. 
Figure O shows p~^dQ/dt as a function of ft for the three slowest sweep rates 



from Fig. |T3|, along with the quasi-equilibrium CV from Fig. y. For the faster 
sweep rates, we used L = 256, while for the slowest sweep rates, L = 128 
was used so that simulations could be performed in a reasonable amount of 
time (about 30 days of CPU time per cycle for p = 3 x 10~^ meV/MCSS on 
a Pentium III 500 MHz processor). To reduce the effects of statistical noise, 
the derivatives were calculated by the Savitzky-Golay smoothing method us- 
ing a quadratic fit to the 0(t) data over a p range of about 60 meV |]58| . 
This accounts for the reduced peak height for p = 0, compared to Fig. ^ 

As expected, the CVs approach the quasi-equilibrium result as p ap- 
proaches zero, and the phase lag between the response and the potential 
depends on the sweep rate. For non-zero sweep rates, there is asymme- 
try between the positive-going (ordering) and negative-going (disordering) 
scans caused by the differences between the phase-ordering and disordering 
processes discussed in Sec. |[ We denote the position of the sharp peak 
in dQ/dt for the positive-going scan as Pp and for the negative-going scan 
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as Jin- Figure. |T3|(a) shows /ip — fic and |/2n — fid as functions of p. The 
peak can be distinguished in the positive-going scans for all of the simulated 
sweep rates but is only visible in negative-going scans for the three slowest 
rates. We therefore report the separation between the positive-going and 
negative-going peak positions, /ip — /in, for only these slowest rates, shown 
in Fig. |15|(b). No simple behavior is seen, but by examining experimental 
peak separations 0, |^, |T0| , it is obvious that our slowest sweep rate is much 
faster than typical experimental sweep rates in the range of 1-10 mV/s. At 
the slowest simulated sweep rate, our peak separation is about 60 mV (in 
experimental electrode-potential units). From the shape of the curve it is 
tempting to suggest that by decreasing the sweep rate by another order of 
magnitude, we might obtain peak separations in the experimental range of 
about 10 mV. 

7 Summary and Conclusion 

We have studied the static and dynamic behavior of Br electrosorption on 
the Ag(lOO) single- crystal surface by Monte Carlo simulation. We have made 
comparison with equilibrium experiments to determine the strength of the 
lateral dipole-dipole repulsion and the electrosorption valency. The simula- 
tions reproduce the experimentally observed second-order phase transition 
between a disordered adlayer and a doubly degenerate c(2 x 2) ordered phase 
at room temperature. As was also found previously by Koper [§, the simu- 
lations reproduce the broad shoulder seen in experimental cyclic voltammo- 
grams. We explain the shoulder by small configurational fluctuations in the 
disordered phase that locally resemble phases with B = 1/4, one of which 
is present for this model in equilibrium at very low temperatures. Stronger 
effects of this low-temperature equilibrium phase might be observed in the 
behavior of other halides on Ag(lOO). 

We also studied the dynamics of phase ordering and disordering pro- 
cesses by stepping the potential between the disordered and ordered phases. 
The ordering process is well described by t^^^ dynamical scaling, consistent 
with growth of domains of the two degenerate c(2 x 2) ordered phases. For 
disordering, no such competition between degenerate domains exists, and 
the process is well described by random desorption with diffusion. Study 
of such disordering processes may provide a way to estimate the adsorp- 
tion/desorption barrier and obtain a relation between MC time and physical 
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time. 

In simulations of CVs for a wide range of potential sweep rates we found a 
pronounced asymmetry between the positive-going and negative-going scans. 
This asymmetry is caused by the differences between the ordering and dis- 
ordering kinetics. No simple behavior is seen for the peak separation as a 
function of the sweep rate, but comparison with experiments suggests that 
our slowest sweep rate is considerably faster than experimental sweep rates 
which are in the range of 1-10 mV/s. 

We leave comparison of simulated and experimental dynamics for future 
work when it hopefully becomes possible to experimentally study the dynam- 
ics of this system in full detail. 
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Appendix 

In this Appendix we discuss the nonlinear fitting of simulated adsorption 
isotherms to the experimental data from Ref. 

Fixing a value of 0nnn, we generated a room-temperature isotherm by MC 
simulation, in this Appendix denoted by 0™(/i;0nnn)- We compared this 
isotherm to three experimental coverage isotherms [Q^^^{E; Ci), Q^^^{E] C2), 
and Q'^^^{E] Cs)] at three different Br~ concentrations [Ci, C2, and C3]. For 
convenience, we set the reference concentration Cq = 1 mM, corresponding 
to the units used for the experimental KBr concentrations. Equation (2) was 
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inverted to obtain E{fl; fiQ,'y,C), assuming trial values of the parameters 
/io and 7, and the /i-axis was rescaled into an E-axis for each concentra- 
tion. Thus, we obtained three simulated approximations to the individual 
experimental isotherms [6'^™(i?; (pnnn, 7, fio', Ci), Q^^^{E; 0nnn, 7, fio] C2), and 
e^^{E; (/.n„n,7,/^o;C3)]. 

In order to obtain the "best" values of the parameters ^nnn, 1, and p,Q, we 
defined the simultaneous squared difference between the simulation at 0nnn 
and the experimental isotherms, 

3 ^max(^') 

X'(0n„n,7,/^O) = E E [0''"(^f;/^O,7,<^nnn;Cfc)-0'=^P(Ef;C,)]^ (9) 
k=l 1=1 

Here Ck is the fc*^ experimental concentration value, Ej' is the value of 
the electrode potential for the k^^ experimental isotherm, J2i is the sum 
over all experimental data points in an isotherm, and X)fc is the sum over 
the three experimental concentrations. The calculation requires that the 
simulation and experimental isotherms be evaluated at the same E values. 
This was accomplished by linear interpolation of the simulation data. 

The parameters 0nnn, 7, and /iq can, in principle, be found by minimiz- 
ing with respect to all three parameters. However, generating a simu- 
lated isotherm for a particular value of 0nnn is too computationally intensive 
to allow simultaneous minimization of with respect to all three param- 
eters. We therefore minimized in two steps, which saves time and al- 
lows better control over the fitting process. First, we generated a family 
of simulation isotherms for 0nnn between —98 and —10 meV at intervals of 
2 meV with a p, resolution of 50 meV over the electrochemical potential range 
—200 to +1000 meV. For each of these simulated isotherms, we minimized 

with respect to 7 and /iq and denoted the minimum error obtained by 
Xmm(0nnn,7™°, Ato"™), whcrc 7™"'(0nnn) and /i™"'(0nnn) are the values of 7 
and /io which minimize for a particular value of 0nnn- In the second step, 
we minimized Xminl'/'nnn, 7™", Z^™™) with respect to (pnnn- The final parame- 
ter results are 0nnn = —26 ± 2 meV, 7 = —0.73 ± 0.03, and /io = 0.67 ± 0.03. 
These errors are estimated from the (pnnn resolution and the differences in 
best-fit parameter values around (pnnn = ~26 meV. The results of the fit are 
shown in Fig. 
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E / V vs SCE 



Figure 1: A single Monte Carlo (MC) simulation isotherm {L = 32, T 
290 K) fit simultaneously to clironocoulometric coverage isotherms at three 
different Br~ concentrations [Q. The parameters, 0nnn = —26 ± 2 meV, 
7 = —0.73 ± 0.03, and fio = 0.67 ± 0.03 meV, were obtained by a non-linear 
fitting procedure described in the Appendix. The reference concentration 
is set to Co = 1 mM, and the electrode potential is given relative to the 
Saturated Calomel Electrode (SCE). 
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Figure 2: MC isotherms for L = 32 and L = 128, T ^ 290 K, and 0nnn = 
—26 meV. A second-order phase transition between a low-coverage disordered 
phase and a c(2 x 2) phase with = 1/2 is observed at /i,. ^ 180 meV. 
©s is the order parameter for the c(2 x 2) phase. In an infinite system, 
I ©si oc {ft — ficY^^ for fj, > pLc- The L = 32 isotherms arc generated from 
10,000 independent samples for each value of Ji. For L = 128, the number 
of independent samples per p, value varies between 1000 (away from the 
transition) and 3000 (near the transition). Thus, the statistical errors are 
larger than for L = 32. 
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Figure 3: Quasi-equilibrium CV for L = 32 and L = 128 at T ^ 290 K 
with 0nnn = — 26 HieV. The derivative dQ/dji is calculated using a three- 
point difference method. The peak corresponds to the second-order phase 
transition between the disordered and c(2 x 2) phases. The broad shoulder 
on the negative side of the sharp peak is due to configurational fluctuations 
of the Br adlayer as discussed in the text. The statistical errors are larger 
for L = 128 than for L = 32. 
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Figure 4: Approximate phase diagram and observed phases for Br elec- 
trodeposited onto Ag(lOO). (a): Approximate phase diagram obtained with 
a two-site heat-bath MC algorithm for L — 16. The horizontal dashed hne 
indicates room temperature. The disordered phase and both stable ordered 
phases arc labeled. The three observed ordered phases are (along with sym- 
metry equivalents) (b): c(2 x 2) with = 1/2 (stable), (c): p(2 x 2) with 
9 = 1/4 (metastable), (d): c(4 x 2) with 9 = 1/4 (stable). 
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Figure 5: The structure factor, S{k), which is proportional to the surface X- 
ray scattering (SXS) intensity. Equihbrium simulation at room temperature 
for L = 256 and p, = +100 meV, averaged over 1,000 independent configu- 
rations, (a): Three-dimensional view of S{k), with the grayscale indicating 
logiQ[S'(fc)]. The delta-function peak in the far corner is the Bragg peak at 
(0,0), corresponding to 6 ~ 1/4. (b): Contour plot of S{k) with the same 
grayscale as in (a). The symbols show where Bragg peaks would occur for 
fully ordered phases [c(2 x 2), x; p(2 x 2), Q; c(4 x 2), □]. The peaks on 
the zone boundary {kx = 2tt, ky = 2tt) are also shown for clarity, (c): A 
32 X 32 portion of a real-space configuration. Small domains resembling all 
three ordered phases shown in Fig. ^ are noticeable. 
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Figure 6: Schematic energy landscape within the Butler- Volmer approxi- 
mation. The states / and F denote lattice-gas states. The state T\ denotes 
the intermediate transition state for a microscopic process of type A. The 
"bare" barrier associated with A is denoted A;^. 
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Figure 7: Time series for sudden disorder-to-order potential-steps at room 
temperature for L = 256 and 0nnn = —26 meV. The insets show the correla- 
tion length for the c(2 x 2) order-parameter fluctuations vs. t^^"^. (a): Deep 
step from /2i = —200 meV to fl2 — +600 meV, averaged over 10 independent 
runs, (b): Shallow step from p, — —200 meV to p,2 — +250 meV, averaged 
over 7 independent runs. 
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Figure 8: Sequence of structure factors and real-space configurations fol- 
lowing a deep disorder-to-order stepv^j^/ii = —200 meV, /Z2 = -1-600 meV, 
L=256, T ?a290 K. (a,c,c): Structure factor in three-dimensional cut-away 
and full projected views at times 25, 2000, and 10,000 MCSS, respectively, av- 
eraged over 10 independent runs. The grayscale on each is the same. (b,d,f): 
128x128 section of a 256x256 real-space configuration at times 25, 2000, and 
10,000 MCSS, respectively. The two different shades of gray idicate Br on 
the two sublattices. 




Figure 9: Circularly averaged structure factor at various times for a deep 
disorder-to-order step from fii = —200 meV to /i2 = +600 meV for L = 256 
and 0nnn = ~26 meV. For comparison with Porod's Law p^, EM, the heavy 
line indicates the power law |fc — /c^ttI"^- 
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Figure 10: Time series for sudden order-to-disorder potential steps at room 
temperature for L — 256 and 0nnn = —26 meV. (a): Deep step from fti — 
+600 meV to p,^ — —200 meV, averaged over 3 independent runs, (b): 
Shallow step from fti = +600 meV to /i2 = +100 meV, averaged over 4 
independent runs. The inset shows the correlation length for the c(2 x 2) 
order-parameter fluctuations vs. t. All four dynamical regimes discussed in 
the text are labeled with Roman numerals. 
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Figure 11: Time series for shallow negative potential step with various 
values of the adsorption/desorption barrier Aa. All other parameters have 
the same values as in Fig. |10|. Time is given in dimensionless units as 
t' = tuexp {—Aa/kBT). The data show the results of a single potential-step 
simulation for each value of Aa- 
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Figure 12: Parametric plot of |0s| vs. for a shallow negative potential-step 
with various values of the adsorption/desorption barrier Aa- The data shown 
are the same as in Fig. |TT|. The arrow indicates the direction of increasing 
time. 
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Figure 13: G vs. /i for various sweep rates p. Dynamic MC simulations 
were used to generate the curves for p 7^ 0. For p=0, the equihbrium cov- 
erage from Fig. || is shown. L=256 for p=lxlO"^ meV/MCSS to 3xl0~^ 
meV/MCSS, L=128 for p=lxlO-=^ and 3x10"^ meV/MCSS, and L=32 for 
p=0. Here, /ii=-200 meV for p=l x IQ-^ meV/MCSS to 1 x lO'^ meV/MCSS 
and pi=— 400 meV otherwise, while p2 = +600 meV. The arrows indicate 
the direction of increasing time. 
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Figure 14: Dynamically simulated cyclic voltammograms (CVs) for various 
sweep rates showing dO/dt vs. fl. For all non-zero sweep rates, the dynamic 
response dQ/dt, which is proportional to the voltametric current, is normal- 
ized by the sweep rate. For p=0, the quasi-equilibrium CV, d9/d/i from 
Fig. ^ is shown. All parameters and symbols are as in Fig. |TB|. 
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Figure 15: Peak positions and peak separations for various sweep rates, 
fa): Positive-going and negative-going peak positions vs. p. All parameters 

Jin- For 
-0.73. 



are the same as in Figs. |T3| and |TJ. (b): Peak separations fip 
comparison with experiments, the units are given in mV using 7 
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Table 1: Initial adsorption, desorption, and nearest and next- nearest neigh- 
bor diffusion rates, Ra, Rd, Run, and Rnnn, respectively, for each of the 
four potential steps. The rates are given in (MCSS)~^. The parameters 
are /cbT = 25 meV, 0nnn = -26 meV, = 300 meV, A„n = 100 meV, 
Annn = 200 meV, and u — 1. 
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